Thermodynamics of the Bose-System with a Small Number of Particles 



O 

o 

(N 



I. 0. Vakarchuk*, A. A. Rovenchak* 
Department for Theoretical Physics of Lviv National University, 
12 Draghomanov St., UA-79005, Lviv, Ukraine 

A theoretical description of the interacting Bose-system is proposed. It is based on the extrap- 
olation of the results obtained for the systems with a small number of particles N = 2, 3, 4, etc. 
to the bulk case of N = oo. It is shown that already the system with N = 12, 13 behaves almost 
as a bulk in a wide temperature range. Special attention is paid to the phase transition in these 
systems. The hard sphere potential is used in calculations. The sequence of heat capacity maxima 
is approximated as C™ ax /iV ~ 13.6 — aN~ £ with e = 0.0608 giving the value of bulk heat capacity 
as 13.6 while experimental value is close to 16. The temperature of A-transition is estimated as 2.1- 
| 2.3 K versus experimental 2.17 K. Quite good qualitative and satisfactory quantitative agreement 

with the experimental data has been achieved. 
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& 1 !• INTRODUCTION 

Liquid 4 He is the only atomic system having its properties determined mainly by quantum effects and the most 
essential among the latter is the statistical effect of identical particles. Due to it the macroscopic number of the 4 He 
atoms can have a zero momentum at some conditions, i. e., it can occupy a single quantum state. The phenomenon 
of bosons accumulation at the lowest energy level is called Bose- Einstein condensation (1924). According to London's 
£h ' assumption (1938), it is believed that superfluidity, A-transition, and other unusual properties of liquid He are a 
consequence of this very effect. 

However, despite the lapse of time that has passed since the discovery the phenomena of A-transition by Keesom and 
Clausius in 1932 and the superfluidity by Kapitsa in 1938 until now noone has managed to create a comprehensive 
-H , microscopic theory explaining the properties of strongly interacting Bose-systems such as 4 He within the whole 
temperature range of < T < oo. 

Bogoliubov |]J was the first to show how one can build a theory of real helium proceeding from the approximate 
second quantization techniques for the model of weakly-imperfect Bose-gas. This theory works well for T < 1 K, 
especially while T — > 0. 

The subsequent theoretical studies brought about the appearance of the quantitative theory of the ground state 
of a strongly-imperfect many-boson system leading to a good agreement with the experiment and the creation of 
the theoretical description of low-energy excitations which are responsible for the low-temperature behaviour of 
thermodynamic functions. The experimental investigations of the spatial structure by means of X-ray and neutron 
diffraction and of the structure of the energy spectrum and the condensate fraction by means of the inelastic neutron 
scattering revealed a fairly good theoretical explanation. 
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1 , I As a consequence of the experimental evidences it was supposed for a long time that He heat capacity had a 
logarithmic divergency at the transition temperature T\. Brout showed in his early works based on the Hartree- 
Fock method that it was impossible to find an explanation of the observed heat capacity behaviour within simple 
theoretical approaches. There is an interesting notice regarding this question Feynman made in his book ||: "The 
explanation of this behavior is left as an exercise for the reader. If successful, publish!" 

• i-h , Feynman suggested that the A-transition should be studied using the full A^-particle density matrix for the ideal 
Bose-gas with a phenomenologically introduced effective mass for particles Q. This approach leads only to a shift 
of the Bose-condensation temperature in comparison with the ideal case leaving the form of the heat capacity curve 

■ - - i unchanged. 

In the early 1970s the studies of the A-transition were started in the renormalization group (RG) approach for 
the theory with the degenerate two-component order parameter. The rigorous grounds of this approach based on the 
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functional integration techniques are given in Refs. g, They use the coherent states representation with the complex 
order parameter the absolute value of which equals -jNo where No is the number of atoms in Bose-condensate. 

In this approach the heat capacity critical exponent was calculated as a — 1/8. It obviously does not correspond 
to the logarithmic divergence for which a equals zero. However, taking into account that the calculated value of a 
is a small magnitude, the experimental confirmation of this result required very precise measurements which became 
available only recently 0. Before it, there was even no certainty about the sign of a. The measurements showed that 
the critical exponent a is negative, i. e., the heat capacity is finite at the transition temperature. 

Recently the theoretical confirmation of this result appeared. Kleinert |||| obtained for a the respective values of 
a = -0.01294 ±0.00060 and a = -0.0120 ± 0.0009, while the experiment gives a = -0.01285 ± 0.00038 §. One can 
find a bit different approach to the critical exponents calculation in the paper by Campostrini et al JTo) ] . The authors 
obtained the value of a — —0.0150 ± 0.0017. But, despite a very good agreement of these results with the experiment, 
RG works only in the immediate vicinity of the transition point T\ and gives no possibility of getting some information 
of the thermodynamic functions behaviour in a wider near-critical range, e. g., for |T — T\ \ ~ 10~ 3 K. 

In addition to these theoretical investigations, the numerical studies of 4 He thermodynamic properties were con- 
ducted. The Diffusion Monte Carlo (DMC) techniques are applied for the ground state calculations, and the Path 
Integral Monte Carlo (PIMC) is used for finite temperatures. These techniques have been applied since the 1960s. 
At present, PIMC is the most effective method. Its only essential imperfection lies in the fact that it is impossible 
to use it for the ground state (T = K). But one can avoid this problem by considering the limit of the results at 
T — > K. A detailed analysis of PIMC might be found in the work by Ceperley jy]]. It appears that the numerical 
approach provides a good agreement with the experimental data in a wide temperature range on the whole, except 
for the vicinity of the A-transition point. 

Thus, we have the situation when no satisfactory theoretical description of 4 He exists for a considerable temperature 
range of 1 K < T < T\ and T > T\, and the only way out is the application of numerical techniques. It gives some 
space for the investigations in this domain. Not claiming the solution of the problem, we still propose here the model 
which, in our opinion, can near us to the understanding of the processes in many-boson system via the consideration 
of systems with a small number of particles. 

For the latter one can expect a much more precise description in terms of statistical mechanics in comparison with 
what might be reached while directly considering a many-particle system. 

The aim of the present work consists in obtaining a thermodynamic description of an interacting many-boson system 
basing on the results for the systems with a small number of particles N. We start our analysis from N — 2, 3, 4, 
etc. in order to extrapolate these results for the case of N = oo after establishing some dependencies on N and the 
temperature. We pay a special attention to the fact that even the system consisting of a small number of particles 
exhibits the peculiarities on the heat capacity curve such as the maximum which means the phase transition in the 
limit of N — > oo. We will show at the same time that the system with N = 12, 13 tends to reveal qualitatively an 
"almost bulk" behaviour in a wide temperature range. 

The boson systems with a relatively small number of particles have been experimentally obtained only in the last 
couple of years as particles in traps (see Refs. 12, 13 and references therein). It gives us a possibility of developing 
some methods for such systems and a further comparison of the results with the experimental data. 

We use the partition function formalism to calculate thermodynamic characteristics of the free Bose-system. A 
method somewhat similar to this one was used in works jl2|,[l3| for the study of the trapped finite Bose-systems: the 
recursion formulae for the partition function proposed by P. Borrmann and G. Frankc 141. But our expressions do 
not take the restriction for the system energy to be written as a sum of one-particle energies. So, one can use them 
not only for the ideal gas but for the interacting system as well. 

The paper is organized as follows. In Section II some known initiating expressions for the calculating the thermo- 
dynamic function are given in short. In Sections III and IV we consider how the proposed method works in the case 
of an ideal Bose-gas the results for which are well-known in the thermodynamic limit when the number of particles 
N — > oo, The next step is the determination of the asymptotic behaviour of our results depending on the number 
of particles N. We will try to expand the found regularities for the case of interacting bosons. The corresponding 
calculations for the hard-sphere potential are given in Section V. 

Generally, it is impossible to obtain a non-analyticity point of the heat capacity function having finite number of 
particles. Therefore, no "pure" phase transition can be detected in this model. But we expect a certain qualitative 
agreement with the experiment to be quite sufficient for the approximation used. 
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II. GENERAL STATEMENTS 



We consider the system of N particles with the Hamiltonian Hn- The partition function Zn — Spe~@ HN where (3 
is inverse temperature, [3 = 1/T. 

In the coordinate representation we have 



Zn = J dr x ... J dr N R N (ri, . . . ,r N | r x , . . .,rjv), 



(1) 



where Rn is the density matrix. The particles are bosons of the mass rn ( 4 He atoms) with their coordinates 
restricted in volume V. In this case Rn is given by [15l 



(2) 




where Q is counting permutations of (1, . . . , TV) and A is the thermal de Broglie wavelength 

-2\ V2 



A 
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The factor Pn takes into account the interaction and will be described further. 
The free energy Fn is given by 



total energy 



and heat capacity 



En — ~qp^Fn), 



c N = p 2 



7" 

N 



z 



N 



Z'n 
Zn 



the primes mean derivation over (3. 



(3) 



(4) 



(5) 



(6) 



III. PARTITION FUNCTION FOR THE SYSTEM OF FREE PARTICLES 



The ideal Bose-gas is a well-studied system for which the expressions for the thermodynamic functions are known. 
Thus, it can be a good test example for the applicability of our method and we are going to use this fact here. 

On calculating the thermodynamic functions of the ideal Bose-system with a small number of particles and com- 
paring them with the well-known results for infinite system one can detect some regularities showing the point when 
several particles start to demonstrate the statistical behaviour. 

One can find an expression for the partition function of N ideal bosons in the explicit form. For this purpose we 
write the density matrix Rn as follows: 



Rn(ti, ...,t n \ti,...,t n ) 



where 



1 1 

m a 5 ^ 



A 



N ■ 



(7) 



1 K\% K13 . . . Kin 
K 2 \ 1 i^23 ■ ■ ■ K 2 n 

Kni K N 2 K N 3 ■ ■ ■ 1 



(8) 



3 



and 

K ij =K{r ij ) = e-^y x \ r« = |r y |, r,, r, v,. (9) 

The subscript "+" means that every "— " in the expression for the determinant should be substituted with "+" . 
The partition function for free particles is 

Zn = Ma^/* 1 ■■■J drN An - ( 10 ) 

The "determinant" (||) might be written as |l5| 

A N = A N ) _ 1 + K Xi K n \%'%+ J2 {K li K ij K jl + K lj K ji K a )A^ ) +..., (11) 

Ki<N Ki<j<N 

The superscripts in An-i mean the lack of the correspondent rows and columns in the expression for Ajy— 
Integrating (|l(]| ) we obtain 



ZN -mi- 



■ Zn-i + j^yZ N - 2 J dvi J dr 2 K 12 K 21 + -^^Z N - 3 J dv x j dv 2 J dr 3 K 12 K 23 K 31 + .. .j , (12) 
After the Fourier transformation 

K <r) = ^ y E e<VK «> ( 13 ) 

q 

Kq = J dr e-^ r K{r) = X^e^^l^ (14) 

the partition function is 

i N Z N _i 

ZN = ^l>2iw> z ° = h (15) 
h i=i 

One can use this expression to calculate the partition function of the system. We will make such calculations in 
the next section. 



IV. MAIN RESULTS FOR THE FREE PARTICLES 



It is possible to calculate heat capacities Cn for N = 1, 2, 3, . . . particles sequentially using expressions (||) and 
(|l5|). In addition, we can show explicitly the low- and high-temperature asymptotics for Cn, namely: 



(16) 
(17) 





3 


If 


T— oo ^ 


Cn 


3 1 


If 





= -— + A N T 3/2 , A N = const with respect to T. 



In the case of N = oo this temperature dependencies coincide with the known results for the ideal Bose-gas [[l6) . 
For low temperatures one has 



N 



15C(5/2) fT \ 



N= 



4 C(3/2) \TJ 



3/2 



(18) 



the high-temperature limit being 3/2. The critical temperature T c (the non-analyticity point on the heat capacity 
curve) in this case is defined from the equation 



p\l = C(3/2), 



(19) 
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where A c is the thermal de Broglie wavelength for T c , and £ is the Rieman ^-function, £(3/2) = 2.612 .... Taking 
p = 0.02185 rA~ 3 for liquid 4 He we obtain T c = 3.138 K. 

In this work we consider the systems with different N from 2 to 1200. As was mentioned above it is impossible 
to obtain the non-analiticity in the heat capacity, so we take the point of Cat maximum as the critical temperature 

It appears that those values could be described by the linear dependence 



P 



AW 
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2 - 61 -^T/3< C ^ L9 ( 2 °) 



with \i N ^ being thermal de Broglie wavelength for Tc N \ One can easily notice that 2.61 is very close to £(3/2) = 
2.612... . 

We also have a similiar correlation for the Cjv maxima C™ ax /7V. So, one should expect that the limit N — > oo gives 
the correct curve (see 





FIG. 2. Temperature of the Cn/N maxima in the ideal case, 
Solid line shows dependence (pp|), 
circles correspond to the data from Table I. 
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The results are also given in Table I. 



AT 


n 

P A c 


/^fmax / AT 

C N /iV 


AT 


P A c 


c N /J\ 


l 




1.500 


90 


2.185 


1.747 


2 


1.131 


1.532 


100 


2.200 


1.753 


3 


1.297 


1.553 


200 


2.286 


1.783 


5 


1.490 


1.582 


300 


2.328 


1.799 


7 


1.606 


1.602 


400 


2.354 


1.809 


10 


1.717 


1.624 


500 


2.373 


1.816 


15 


1.830 


1.649 


600 


2.387 


1.822 


20 


1.902 


1.666 


700 


2.399 


1.827 


30 


1.993 


1.690 


800 


2.408 


1.831 


40 


2.050 


1.706 


900 


2.416 


1.834 


50 


2.091 


1.718 


1000 


2.423 


1.837 


60 


2.123 


1.727 


1100 


2.429 


1.839 


70 


2.147 


1.735 


1200 


2.434 


1.842 


80 


2.168 


1.742 


oo 


2.612 


1.926 



TABLE I. Main results for the ideal systems 



As one can see from Fig. 1, good qualitative and satisfactory quantitative agreement is reached at N = 100. 
Although, typical elements of the heat capacity curve appear at smaller N, e. g., a clearly defined maximum is found 
at N — 10-15. As we shall show further, these very values might be considered as the limit between "individual" and 
statistic behaviour of the system with N particles. 



V. INTERACTION IN THE HARD-SPHERES POTENTIAL APPROXIMATION 



In a many-body system, it is a very complicated problem to take into account the interaction and it has not as 
yet been completely solved. Specifically, the systems with a small number of particles give a possibility to consider 
this problem from the first principles, i. e., to include the interaction into the calculations explicitly. From the 
mathematical point of view, the computations become heavily complicated while the number of particles increases, 
and the direct interaction account soon turns out to be impossible. 

We will proceed from the expression for the density matrix in the free particles case. In the systems like 4 He 
small interatomic distances are essential when atoms behave almost as hard spheres. Thus, the wavefunction should 
vanish rapidly as soon as the distance between any two atoms approaches the hard sphere diameter. 

Due to this we shall take the interaction into account by the substitution: of An in (|l(]) with AnPn where Pn is 
a function of the particles coordinates: 



dn . . . / dr N A N P N . (21) 



7 _ 1 1 
Zn ~ A!A^v 

For the sake of simplicity we suppose the interaction to have a pair character since the obtained results are expected 
to be rather qualitative than quantitative. In the case of hard spheres potential we have 

JV-l N 

P N (r u ...,r N ) = [] l[P(i,j) 7 

i=l j=2 

P(i,j)=P(\r i -r j \)=P(r ij ), 

p(i,j) = i+/(i,i)=| 5; Zi l ' (22) 



a is the hard sphere diameter. 

One can write the partition function as 
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1 N 

Z N = ^Y J b' l {N)Z N _ l , Z Q 



(23) 



where 



is the effective mass of He atom, 



2irf3h' 



1.7m ill 



1/2 



(24) 



We integrate ( ^lj ) in two steps. The first one is the circular approximation ( CA ) and the second one is the free-volume 
approximation (FVA). 

1) CA is to the effect that only the interparticle interactions in the circle 1 — 2 — 3 — . . . — N — 1 are taken into 
account. It means we neglect the interaction between the first and the third, fourth, etc. particles, between the 
second and fourth, fifth, etc. particles, .... This makes it possible to integrate ( |2~l] ) using Fourier transformation as 
we did in (JTc|) . Using also FVA we obtain for the coefficients 6J(JV): 



b[(N) = < 

Here we have used such designations: 



[l-f(N-l)Y[l-f(N-3)} II (iV -*)]/,_!, ;>3 

k=3 



[1 



I < 3 



(25) 



(26) 



In = < 



8(4tt) 



n-l 



2 | _L p -QV47r _ 
47T 

1. 



a/A 



x sin (Jx e 



\ n+l 
" rx2 dx 1 , n > 



(27) 



n = 



2) In FVA we indirectly take into account the interactions neglected in CA. For this purpose each integral over r», 
i = is multiplied by the factor (1 — nfo/U) where n equals to the number of multipliers P(i,j) which are 

not taken into account in the Fourier transformation. To clarify the above statement we consider the most simple 
non-trivial example with N = 4. According to (E5I) we have: 



&i(4) = 



i - 



i 



3-^q 
V 



Iq. Integral over ri "hooks" the folloving three multipliers: P(l,2), P(l,3), P(l,4). 



2vo 
V 



h. Integral over n "hooks" P(l, 2), P(l, 3), P(l, 4) but P(l, 2) is included into h. Integral 
over Y2 "hooks" more P(2,3) and P(2,4). 

3 

12- Each integral over "hooks" additional P(i, 2). 



1 



V 



1 



V 



I3. P(l,3) and P(2,4) are not included into I3. 



Thus, FVA improves the circular approximation. But it is not possible to use it for large / and N because in this 
case a big number of P(i,j) multipliers will not be integrated in the proper way but via formal factors. 

In the calculations we assume the hard sphere diameter a — 2.1 rA. The heat capacity curves exhibit a fast 
approaching of the maximum point to a value from the region of 2.2-2.3 while N is increasing. The exact A-transition 
temperature is 2.17 K. 

As has already been mentioned, in the ideal case the linear (over N~ e ) dependence fits the sequence of Cn/N 
maxima quite well 



f max 

N 



CT X - aN- 



e = 1/3. 



(28) 
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Assuming the same form for the hard spheres system one can find C™ ax = f 3.6 and e = 0.0608 making the maxima 
sequence as follows: 



N 



13.6- 12.9 N~ £ , £ = 0.0608. 



(29) 



Although small N values of 12, 13 do not allow one to ide ntify the behaviour at N — > oo we obtain the values of 
heat capacity maximum from the range of experimental data EaJlj^ M . 



~i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 r 



c 



N 



N 



— N = 


18 


— N = 


13 


-N = 


G 


- N = 


2 



NIST 




1 2 3 4 5 

FIG. 3. Heat capacities in the free volume and circular approximation. The circles show the NIST Cy data ]lg| ]. 



One can write the exact expression for the partition function similiar to that in the case of CA and FVA: 

N 



N — l: 



Z = l. 



1=1 



From expansion (|Tl]) and formula (|2l]) we obtain (without approximations): 

1 1 1 



bKN) 



V(N-l)\ (A 3 )^- 1 
i-i i 

x ii n p (^') b IN (i,...,i) 

i=l j=i-\-l 



dr 1 ... / dri K 12 . ..Kt-ijKn 



(30) 



(31) 



B lN (l, . . . ,1) = / dri+i 



A (i,...,0 p (i,...,0 



N 



dr 



N- 



n n p{ij) 
i=ij=i+i 



Zn-i 



(32) 



(i,...,i) 

N-l 



JV-1 N 

n n p<m)- 



(33) 



We also introduce the following designations to simplify further records: 



R(l,2) = / dr 3 P(1,3)P(2,3), 



K ij = K ij P{i,j). 

One can write the coefficients tf(N) independently of the interatomic potential, i. c. for any Pn 

*i(i) = i. 

6 i( 2 ) = ^2 / dr i / *2 P(l,2), 



61(3) = 



dr 2 P(1,2) J R(1,2) (l + tf 12 tf 21 ) 



y rfn /" dr 2 P(l,2) (l + K 12 K 21 ) 

65(2) = / dl "2 ^12^21^(1,2), 

6 2(3) = y dri y dr 2 K 12 K 21 P(1, 2)R(1, 2), 

6 a(3) = "^^2 y dri y ^r 2 y dr 3 K 12 K 23 K 31 . 

These expressions might be written via one-dimensional integrals: 

Air r°° 

K(2) = — J P(r 12 )r 2 12 dr 12 , 



47T 

A3 



/ K 2 12 P(r 12 )r 2 12 dr 12 , 
Jo 



8 f°° dq ( f°° 1 

R(ri 2 ) = — / — singri 2 < / rP(r) sin qr dr > , 
r i2 J q [Jo J 

63( 3 ) = / — (/ ^i2-P(fi2)r-i 2 sinqr 12 dr 12 | , 



(A 3 ) 2 7o 9 Uo 
65(3) =yjjij K 2 12 P(r 12 )R(n 2 )r 2 12 dr 12 . 

In the hard spheres approximation integration is made simply. We obtain 

6*(1) = 1, 
«(2)«l-£, 

«(2)-^{l-f(^} + K^, 

= V-3v + §i + X%(3) 
l[ ' Vb\{2) + \ 3 b* 2 (2) 



12V i V. A ; 7r A 2 27r 2 y ' \ A 4 7r A •_>.-- 

^{i(- rf (^)M 



65(3) = 32^y°° ^ |^: e_Q2/4 " - y^^ sin ^ ^ dx 



Our R coincides with the correspondent function from p0| : 

R(ri 2 ) = J dr 3 (1 + / 13 )(1 + / 23 ) - V - 2v + 9l (r l2 ), 
The function gi(r\2) is given by 

( 



51 (r) = / dr 3 /i 3 / ; 



23 



i^ 3 ( 2 -^ + ^«)"J' r<2a 

0, r > 2a 



(56) 



(57) 



Meanwhile we cannot write exact expressions for greater N. The comparison of the exact results and those in CA 
and FVA is presented in Table II. 



N 


Approximation 


Exact values 


T w 




T w 


C^/N 


1 




1.50 




1.50 


2 


2.86 


1.67 


2.86 


1.67 


3 


2.54 


1.78 


2.19 


1.83 


4 


2.28 


1.90 






5 


2.25 


2.01 






6 


2.24 


2.11 






7 


2.23 


2.20 






8 


2.23 


2.29 






9 


2.23 


2.37 






10 


2.23 


2.44 






11 


2.23 


2.50 






12 


2.23 


2.56 






13* 


2.23 


2.61 






18 


2.23 


2.82 






oo 


2.2-2.3 


13.6** 







We suppose that the system of 12-13 particles behaves almost as a bulk since the atoms in 4 He are packed 
compactly enough to say that one atom is surrounded by 12 atoms. Therefore, 13 atoms constitute the first 
closed sphere (see Fig. 5). 

The experimental value is ~ 16. The heat capacity in the A-point was believed to be infinite for a long time but 
recent experiments show the finite peak M. 



TABLE II. Results for the interacting systems 



14.0 
13.8 
13.6 
13.4 
13.2 

13.0; 

3.0 
2.5 
2.0 
1.5 
1.0 



1 1 1 1 /A 1 1 I 1 


! i i i i 1 i i i i 1 i i i i 








\. o - 


1 , 1 /V , , . , 


. i 



0.00 0.02 0.040.80 



0.85 



0.90 



0.95 



1.00 
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FIG. 4. C^^/N values for the hard spheres system. 
Solid line shows dependence (p^), 
circles correspond to the data from Table II. 




FIG. 5. Helium-4 atoms as hard spheres. 

(a) The first closed sphere of 12 atoms surrounding 1 central atom. One atom is not seen since it is situated exactly behind 
the central atom. 

(b) Solid line shows the "surface of atom" , dashed line shows the "surface of alien atoms" , dashed-dotted line shows "surface 
of inaccessible space". Thus, assuming the volume of "inaccessible space" 4-kR 3 /3 to be the volume occupied by 1 atom one 
can obtain the hard sphere diameter a ~ 2.2 rA. 

Our expressions for the partition function have the correct form zffl in the classical limit (fi — » 0) plj . 

4 ci) =^T3^/^.../^ n p M- (58) 

J J l<i<j<N 

It means that the proposed method is applicable for a wide temperature range 



CONCLUSION 



It is shown that proposed method leads to a good agreement between the calculated results and the experimental 
data. We obtained not only qualitive but even satisfactory quantitive fit as well. 

The expansions over 1 /N obtained in this work demonstrate the existence of some parameter e having the value 
of 1/3 for the ideal system. We estimated e = 0.0608 from the numerical analysis of the results. The theoretical 
calculation of this parameter will be our next problem along this line of research. 
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